Preparing the data

Loading the packages and data directories

options(width=108)
out=F
# The list of subjects, the order of conditions, and the thresholds are derived from Subjects.xlsx
library(xlsx)
library(ggplot2)
library(plyr)
library(reshape)
library(matrixStats)
library(gridExtra)
library(lme4)
library(lmerTest)
library(BayesFactor)
#library(splines)
db <- '/home/egor/Dropbox/' # on Linux
db <- '/Users/Egor/Dropbox/' # Windows
# db <- '~/Dropbox/' # on Mac
source(paste(db, 'Prog/R/myFunctions/pvalfn.R', sep=''))

Plot variables

# theme for plotting:
alpha <- .7
w <- .56 # proportion width in group plots
themefy <- function(p) {
    p <- p + theme_bw() + 
         theme(panel.grid.minor.x=element_blank(), panel.grid.minor.y=element_blank(),
            axis.text=element_text(size=8), axis.title=element_text(size=9),
            legend.text=element_text(size=8), legend.title=element_text(size=9),
            legend.key = element_blank(), legend.margin=margin(t=-.04, unit='in'),
            legend.background = element_rect(fill='transparent'),
            plot.title=element_text(face='bold'))
}
cc <- c('#e31a1c','#fdbf6f','#a6cee3','#1f78b4','#b2df8a','#33a02c','#fb9a99')
xLab <- 'Target Peak Time (s)'
yLab <- 'Log Contrast Threshold'
# colLab <- expression(paste('\nTarget\nVelocity (', degree, '/s)', sep=''))
colLab <- expression(paste('Target Eccentricity (', degree, ')', sep=''))
# colLabType <- 'Mask Type'
yLim <- c(-1.6,-0.6)
dodge <- position_dodge(width=0.1)
sumFn <- function(ss, subjStr='subj', xStr='targTpeak', grpStr='targEcc'){
    sumSubj <- ddply(ss, c(subjStr, xStr, grpStr), summarise,
                     mnS=mean(threshMean), se=sd(threshMean)/sqrt(length(threshMean))) #, .drop=F)
    sumSubjMn <- ddply(ss, c(subjStr), summarise, mnStot=mean(threshMean)) # total mean across conditions per subj
    sumSubj <- merge(sumSubj, sumSubjMn)
    sumSubj$normS <- - sumSubj$mnS / sumSubj$mnStot # normalized subject mean
    sumSubj$seNorm <- NA
    # sumSubj$mnS[is.na(sumSubj$mnS)] <- 0
    sumGrp <- ddply(sumSubj, c(xStr, grpStr), summarise,
                  mn=mean(mnS), se=sd(mnS)/sqrt(length(mnS)),
                  norm=mean(normS), seNorm=sd(normS)/sqrt(length(normS)))
    sumGrp$subj <- 'average'
    sumSubj <- rename(sumSubj, c(mnS='mn',normS='norm'))
    sumComb <- rbind(sumGrp, subset(sumSubj, select=-mnStot))
    sumComb$se[is.na(sumComb$se)] <- 0
    sumComb
}
plotAve <- function(pss, subjStr='subj', xStr='targTpeak', grpStr='targEcc', 
                    xlab=xLab, ylab=yLab, collab=colLab, yStr='mn', seStr='se'){
    pss$yMin <- pss[,yStr] - pss[,seStr]
    pss$yMax <- pss[,yStr] + pss[,seStr]
    pss[,grpStr] <- factor(pss[,grpStr])
    p <- ggplot(pss, aes_string(x=xStr, y=yStr, colour=grpStr, group=grpStr,
                            ymin='yMin', ymax='yMax')) +
        geom_point(position=dodge, size=1, alpha=alpha) + geom_line(position=dodge, alpha=alpha) +
        scale_x_continuous(breaks=c(.5,1,1.5), labels=c('0.5','1','1.5')) +
        geom_linerange(position=dodge, show.legend=F, alpha=alpha) +
        labs(x=xlab, y=ylab, colour=collab) + #ylim(0,1) + 
        guides(colour=guide_legend(keyheight=.3, default.unit='inch'))
    p <- themefy(p)
}
plotIndiv <- function(pss, subjStr='subj', xStr='targTpeak', grpStr='targEcc', 
                    xlab=xLab, ylab=yLab, collab=colLab, yStr='mn', seStr='se'){
    pss$yMin <- pss[,yStr] - pss[,seStr]
    pss$yMax <- pss[,yStr] + pss[,seStr]
    pss[,grpStr] <- factor(pss[,grpStr])
    p <- ggplot(pss, aes_string(x=xStr, y=yStr, colour=grpStr, group=grpStr,
                            ymin='yMin', ymax='yMax')) + 
        facet_wrap( ~ subj, ncol=4) +
        geom_point(position=dodge, size=1, alpha=alpha) + geom_line(position=dodge, alpha=alpha) +
        geom_linerange(position=dodge, show.legend=F, alpha=alpha) +
        scale_x_continuous(breaks=c(.5,1,1.5), labels=c('0.5','1','1.5')) +
        # scale_y_continuous(breaks=c(0,.25,.5,.75,1), labels=c('0','','0.5','','1'), limits=c(0,1)) +
        labs(x=xlab, y=ylab, colour=collab) + 
        guides(colour=guide_legend(keyheight=.3, default.unit='inch'))
    p <- themefy(p)
}

Loading the data

allDataDir <- paste(db,'Projects/mc/data_cfs/cfs',sep='')
dataDirs <- dir(allDataDir)
dataDirs <- dataDirs[grep('cfs',dataDirs)]
colsOfInt <- c('participant', 'dom', 'session', 'mcBv', 'targTpeak', 'targXoff2', 
               'targV', 'stairStart', 'meanRev6')
df <- data.frame()
dfRevs <- data.frame()
dfIntn <- data.frame()
for(curDir in dataDirs){
    print(curDir)
    curDf <- read.csv(paste(allDataDir,'/',curDir,'/',curDir,'.csv',sep=''))
    curDf <- curDf[,colsOfInt]
    curDf$cond <- substr(curDir,19,22)
    subjStairs <- dir(paste(allDataDir,'/',curDir,'/',sep=''))
    subjStairs <- subjStairs[grep('.tsv',subjStairs)]
    for(curStairFN in subjStairs){
        curStair <- paste(allDataDir,'/',curDir,'/',curStairFN,sep='')
        curRevs <- read.table(curStair, skip=1, nrows=1)
        curIntn <- read.table(curStair, skip=4, nrows=2)
        curInfo <- readLines(curStair)
        curInfoCols <- data.frame(subj=curDf$participant[1], dom=curDf$dom[1], 
                                  sess=curDf$session[1], stairStart=curInfo[33], 
                                  mcBv=curInfo[41], targTpeak=curInfo[43],
                                  targEcc=curInfo[37], targV=curInfo[31])
        curDfRevs <- cbind(curInfoCols, curRevs[,2:11])
        nTrials <- ncol(curIntn)-1
        curDfIntn <- curInfoCols[rep(seq_len(nrow(curInfoCols)), each=nTrials),]
        curDfIntn$trial <- seq(1,nTrials)
        curDfIntn$intn <- as.numeric(curIntn[1,2:(nTrials+1)])
        curDfIntn$resp <- as.numeric(curIntn[2,2:(nTrials+1)])
        rownames(curDfIntn) <- NULL
        dfRevs <- rbind(dfRevs, curDfRevs)
        dfIntn <- rbind(dfIntn, curDfIntn)
    }
    df <- rbind(df, curDf)
}
## [1] "mc2_cfs_cent_p12_s2_2017-11-02_1537"
## [1] "mc2_cfs_cent_p13_s2_2017-11-01_1030"
## [1] "mc2_cfs_cent_p17_s1_2017-10-26_1039"
## [1] "mc2_cfs_cent_p18_s1_2017-10-26_1443"
## [1] "mc2_cfs_dyna_p10_s2_2017-10-27_1511"
## [1] "mc2_cfs_dyna_p11_s1_2017-11-02_1439"
## [1] "mc2_cfs_dyna_p14_s1_2017-10-27_1235"
## [1] "mc2_cfs_dyna_p4_s2_2017-10-26_1342"
## [1] "mc2_cfs_peri_p12_s1_2017-10-26_1539"
## [1] "mc2_cfs_peri_p13_s1_2017-10-31_1032"
## [1] "mc2_cfs_peri_p17_s2_2017-11-02_1018"
## [1] "mc2_cfs_peri_p18_s2_2017-11-02_1327"
## [1] "mc2_cfs_stat_p10_s1_2017-10-27_1030"
## [1] "mc2_cfs_stat_p11_s2_2017-11-02_1630"
## [1] "mc2_cfs_stat_p14_s2_2017-10-27_1319"
## [1] "mc2_cfs_stat_p4_s1_2017-10-26_1234"
ds <- rename(df, c(participant='subj', session='sess', meanRev6='thresh',
                   mcBv='maskV', targXoff2='targEcc'))
ds$targEcc <- round(ds$targEcc / 35,1)
ds$targV <- round(ds$targV / 3.5,1)
ds$maskType <- ''
ds$maskType[ds$maskV==0] <- 'static'
ds$maskType[ds$maskV==10] <- 'standard'
ds$maskType[ds$maskV==60] <- 'fast'
# ds$maskV <- round(ds$maskV * 60 / 35, 1)
# ds$maskV[ds$maskV<0.05] <- 0
ds$condSplit <- ''
ds$condSplit[ds$cond=='stat' | ds$cond=='dyna'] <- 'v'
ds$condSplit[ds$cond=='cent' | ds$cond=='peri'] <- 'ecc'
head(ds)
##   subj dom sess maskV targTpeak targEcc targV stairStart    thresh cond maskType condSplit
## 1   12   1    2    60       1.5     0.8   4.6          0 -1.378333 2_20     fast          
## 2   12   1    2     0       1.5     0.8   4.6         -2 -1.375000 2_20   static          
## 3   12   1    2    60       0.5     0.8   0.0         -2 -1.263333 2_20     fast          
## 4   12   1    2    60       1.0     0.8   4.6         -2 -1.341667 2_20     fast          
## 5   12   1    2    10       1.5     0.8   0.0          0 -1.363333 2_20 standard          
## 6   12   1    2     0       1.5     0.8   0.0          0 -1.636667 2_20   static

Transformed data sets

By target

thresh <- ddply(ds, .(subj,dom,sess,maskV,maskType,targTpeak,targEcc,targV,
                      cond,condSplit), summarise, threshMean = mean(thresh),
                      threshSD = sd(thresh))
head(thresh)
##   subj dom sess maskV maskType targTpeak targEcc targV cond condSplit threshMean   threshSD
## 1    4   1    1     0   static       0.5     0.8     0 _201            -1.565833 0.11667262
## 2    4   1    1     0   static       0.5     2.9     0 _201            -1.122500 0.14024284
## 3    4   1    1     0   static       1.0     0.8     0 _201            -1.219167 0.02474874
## 4    4   1    1     0   static       1.0     2.9     0 _201            -0.910000 0.30876996
## 5    4   1    1     0   static       1.5     0.8     0 _201            -1.342500 0.08367430
## 6    4   1    1     0   static       1.5     2.9     0 _201            -1.042500 0.13552880

Plotting across target eccentricities

Static target & static mask

ss <- sumFn(thresh[(thresh$maskType=='static' & thresh$targV==0),])
ssAve <- ss[ss$subj=='average',]
ssIndiv <- ss[ss$subj!='average',]
pTargStatMaskStat <- themefy(plotAve(ssAve)) + ylim(yLim)

Individual plots

(plotIndiv(ssIndiv) + ylim(-2,0))

Stair plots

iss <- dfIntn[dfIntn$mcBv==0 & dfIntn$targV==0,]
(p <- ggplot(iss, aes(x=trial, y=intn, colour=targEcc, linetype=stairStart)) +
    facet_wrap(subj ~ targTpeak, ncol=6) + geom_point() + geom_line() + ylim(-2,0))

Moving target and static mask

ss <- sumFn(thresh[thresh$maskType=='static' & thresh$targV>0,])
ssAve <- ss[ss$subj=='average',]
ssIndiv <- ss[ss$subj!='average',]
pTargDynaMaskStat <- themefy(plotAve(ssAve)) + ylim(yLim)

Individual plots

(plotIndiv(ssIndiv) + ylim(-2,0))

Stair plots

iss <- dfIntn[dfIntn$mcBv==0 & dfIntn$targV==16,]
(p <- ggplot(iss, aes(x=trial, y=intn, colour=targEcc, linetype=stairStart)) +
    facet_wrap(subj ~ targTpeak, ncol=6) + geom_point() + geom_line() + ylim(-2,0))

Static target & standard mask

ss <- sumFn(thresh[thresh$maskType=='standard' & thresh$targV==0,])
ssAve <- ss[ss$subj=='average',]
ssIndiv <- ss[ss$subj!='average',]
pTargStatMaskSlow <- themefy(plotAve(ssAve)) + ylim(yLim)

Individual plots

(plotIndiv(ssIndiv) + ylim(-2,0))

Stair plots

iss <- dfIntn[dfIntn$mcBv==10 & dfIntn$targV==0,]
(p <- ggplot(iss, aes(x=trial, y=intn, colour=targEcc, linetype=stairStart)) +
    facet_wrap(subj ~ targTpeak, ncol=6) + geom_point() + geom_line() + ylim(-2,0))

Dynamic target and standard mask

ss <- sumFn(thresh[thresh$maskType=='standard' & thresh$targV>0,])
ssAve <- ss[ss$subj=='average',]
ssIndiv <- ss[ss$subj!='average',]
pTargDynaMaskSlow <- themefy(plotAve(ssAve)) + ylim(yLim)

Individual plots

(plotIndiv(ssIndiv) + ylim(-2,0))

Stair plots

iss <- dfIntn[dfIntn$mcBv==10 & dfIntn$targV==16,]
(p <- ggplot(iss, aes(x=trial, y=intn, colour=targEcc, linetype=stairStart)) +
    facet_wrap(subj ~ targTpeak, ncol=6) + geom_point() + geom_line() + ylim(-2,0))

Static target & fast mask

ss <- sumFn(thresh[thresh$maskType=='fast' & thresh$targV==0,])
ssAve <- ss[ss$subj=='average',]
ssIndiv <- ss[ss$subj!='average',]
pTargStatMaskFast <- themefy(plotAve(ssAve)) + ylim(yLim)

Individual plots

(plotIndiv(ssIndiv) + ylim(-2,0))

Stair plots

iss <- dfIntn[dfIntn$mcBv==60 & dfIntn$targV==0,]
(p <- ggplot(iss, aes(x=trial, y=intn, colour=targEcc, linetype=stairStart)) +
    facet_wrap(subj ~ targTpeak, ncol=6) + geom_point() + geom_line() + ylim(-2,0))

Dynamic target and fast mask

ss <- sumFn(thresh[thresh$maskType=='fast' & thresh$targV>0,])
ssAve <- ss[ss$subj=='average',]
ssIndiv <- ss[ss$subj!='average',]
pTargDynaMaskFast <- themefy(plotAve(ssAve)) + ylim(yLim)

Individual plots

(plotIndiv(ssIndiv) + ylim(-2,0))

Stair plots

iss <- dfIntn[dfIntn$mcBv==60 & dfIntn$targV==16,]
(p <- ggplot(iss, aes(x=trial, y=intn, colour=targEcc, linetype=stairStart)) +
    facet_wrap(subj ~ targTpeak, ncol=6) + geom_point() + geom_line() + ylim(-2,0))

Group plots

grid.arrange(pTargStatMaskStat + #theme(legend.position='none') + 
                 ggtitle('a: static mask, static target'),
             pTargDynaMaskStat + ggtitle('d: static mask, dynamic target'),
             pTargStatMaskSlow + ggtitle('b: standard mask, static target'),
             pTargDynaMaskSlow + ggtitle('e: standard mask, dynamic target'),
             pTargStatMaskFast + ggtitle('c: fast mask, static target'),
             pTargDynaMaskFast + ggtitle('f: fast mask, dynamic target'),
             ncol=2, nrow=3) #widths=c((1-w)/2,w/2))
## Warning: Removed 2 rows containing missing values (geom_linerange).

Analyses

Scaling function

cent <- function(v){
    v <- apply(v,2,function(x){
        x <- x - mean(unique(x))
        x <- x / max(x)
    })
}

Centered data set

dsc <- ds
centCols <- c('dom','targTpeak','stairStart','targEcc','maskV','targV')
dsc[,centCols] <- cent(dsc[,centCols])
# dsc$maskV <- (dsc$maskV / max(dsc$maskV)) #*2 - 1
dsc$maskV_c <- (dsc$maskV + 1) / 2
dsc$targV_c <- (dsc$targV + 1) / 2
# dsc$targEcc_c <- (dsc$targEcc + 1) / 2
dsc$targEcc_c <- round((dsc$targEcc + 1) / 2,0)
dsc$targTpeak_c <- (dsc$targTpeak + 1) / 2
# dsc$subj <- as.factor(dsc$subj)
dsc$maskType = factor(dsc$maskType, c('static','standard','fast'))
head(dsc)
##   subj dom sess      maskV targTpeak targEcc targV stairStart    thresh cond maskType condSplit   maskV_c
## 1   12   1    2  1.0000000         1      -1     1          1 -1.378333 2_20     fast           1.0000000
## 2   12   1    2 -0.6363636         1      -1     1         -1 -1.375000 2_20   static           0.1818182
## 3   12   1    2  1.0000000        -1      -1    -1         -1 -1.263333 2_20     fast           1.0000000
## 4   12   1    2  1.0000000         0      -1     1         -1 -1.341667 2_20     fast           1.0000000
## 5   12   1    2 -0.3636364         1      -1    -1          1 -1.363333 2_20 standard           0.3181818
## 6   12   1    2 -0.6363636         1      -1    -1          1 -1.636667 2_20   static           0.1818182
##   targV_c targEcc_c targTpeak_c
## 1       1         0         1.0
## 2       1         0         1.0
## 3       0         0         0.0
## 4       1         0         0.5
## 5       0         0         1.0
## 6       0         0         1.0

GLM

pvalfn(lmer(thresh ~ dom + stairStart + sess + maskType * targTpeak_c * targEcc_c +
                 (1|subj), data=dsc[dsc$targV==-1,]))
##                                          estm    se      df    tVal  pVal star
## (Intercept)                            -1.248 0.119  10.910 -10.497 0.000  ***
## dom                                     0.048 0.102   5.933   0.474 0.653     
## stairStart                              0.055 0.011 266.956   5.086 0.000  ***
## sess                                   -0.141 0.031 271.534  -4.631 0.000  ***
## maskTypestandard                        0.457 0.059 266.956   7.691 0.000  ***
## maskTypefast                            0.292 0.059 266.956   4.914 0.000  ***
## targTpeak_c                             0.005 0.065 266.956   0.083 0.934     
## targEcc_c                               0.380 0.059 266.956   6.397 0.000  ***
## maskTypestandard:targTpeak_c           -0.125 0.092 266.956  -1.358 0.176     
## maskTypefast:targTpeak_c               -0.088 0.092 266.956  -0.954 0.341     
## maskTypestandard:targEcc_c             -0.213 0.084 266.956  -2.534 0.012    *
## maskTypefast:targEcc_c                 -0.220 0.084 266.956  -2.612 0.010    *
## targTpeak_c:targEcc_c                  -0.110 0.092 266.956  -1.199 0.232     
## maskTypestandard:targTpeak_c:targEcc_c  0.187 0.130 266.956   1.439 0.151     
## maskTypefast:targTpeak_c:targEcc_c      0.234 0.130 266.956   1.799 0.073    .
pvalfn(lmer(thresh ~ dom + stairStart + sess + maskType * targTpeak_c * targEcc_c *
                 targV_c + (1|subj), data=dsc))
##                                                  estm    se      df    tVal  pVal star
## (Intercept)                                    -1.391 0.104  10.204 -13.378 0.000  ***
## dom                                             0.035 0.091   6.000   0.388 0.712     
## stairStart                                      0.044 0.008 543.000   5.353 0.000  ***
## sess                                           -0.046 0.016 543.000  -2.813 0.005   **
## maskTypestandard                                0.457 0.063 543.000   7.242 0.000  ***
## maskTypefast                                    0.292 0.063 543.000   4.627 0.000  ***
## targTpeak_c                                     0.005 0.069 543.000   0.078 0.938     
## targEcc_c                                       0.380 0.063 543.000   6.024 0.000  ***
## targV_c                                         0.252 0.063 543.000   3.989 0.000  ***
## maskTypestandard:targTpeak_c                   -0.125 0.098 543.000  -1.279 0.202     
## maskTypefast:targTpeak_c                       -0.088 0.098 543.000  -0.899 0.369     
## maskTypestandard:targEcc_c                     -0.213 0.089 543.000  -2.386 0.017    *
## maskTypefast:targEcc_c                         -0.220 0.089 543.000  -2.460 0.014    *
## targTpeak_c:targEcc_c                          -0.110 0.098 543.000  -1.129 0.260     
## maskTypestandard:targV_c                       -0.140 0.089 543.000  -1.566 0.118     
## maskTypefast:targV_c                            0.082 0.089 543.000   0.924 0.356     
## targTpeak_c:targV_c                            -0.017 0.098 543.000  -0.175 0.861     
## targEcc_c:targV_c                              -0.194 0.089 543.000  -2.176 0.030    *
## maskTypestandard:targTpeak_c:targEcc_c          0.187 0.138 543.000   1.355 0.176     
## maskTypefast:targTpeak_c:targEcc_c              0.234 0.138 543.000   1.694 0.091    .
## maskTypestandard:targTpeak_c:targV_c            0.090 0.138 543.000   0.650 0.516     
## maskTypefast:targTpeak_c:targV_c                0.100 0.138 543.000   0.720 0.472     
## maskTypestandard:targEcc_c:targV_c              0.161 0.126 543.000   1.278 0.202     
## maskTypefast:targEcc_c:targV_c                  0.165 0.126 543.000   1.310 0.191     
## targTpeak_c:targEcc_c:targV_c                   0.082 0.138 543.000   0.592 0.554     
## maskTypestandard:targTpeak_c:targEcc_c:targV_c -0.136 0.196 543.000  -0.697 0.486     
## maskTypefast:targTpeak_c:targEcc_c:targV_c     -0.185 0.196 543.000  -0.946 0.345